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Liquid  water  flooding  in  micro  gas  channels  is  an  important  issue  in  the  water  management  of  poly¬ 
mer  electrolyte  fuel  cells  (PEFCs).  However,  in  most  previous  numerical  studies  liquid  water  transport 
in  the  gas  channels  (GC)  has  been  simplified  by  the  mist  flow  assumption.  In  this  work,  we  present 
a  two-phase  flow  model  for  the  cathode  side  of  a  PEFC.  The  GC  is  assumed  to  be  a  structured 
porous  medium  with  the  porosity  of  1.0.  The  two-phase  Darcy’s  law  is  applied  to  both  diffusion  lay¬ 
ers  and  GC.  Based  on  the  developed  model,  the  liquid  water  flooding  in  the  GC  and  its  impact  on 
the  liquid  water  distribution  in  the  diffusion  layers  are  explored  in  detail.  Furthermore,  we  study 
the  effect  of  the  immobile  saturation  on  the  predicted  liquid  water  distribution  in  the  diffusion 
layers.  The  results  show  that  neglecting  the  GC  flooding  leads  to  an  incorrect  prediction  of  liquid  water 
distribution  in  the  diffusion  layers  and  an  overestimation  of  the  cell  performance.  The  gas  flow  rate  in 
the  GC  can  be  optimized  to  achieve  the  best  cell  performance.  Finally,  when  considering  the  immobile 
saturation  in  the  model,  more  liquid  water  is  predicted  in  the  diffusion  layers. 
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1.  Introduction 

In  the  pursuit  of  reduced  dependence  on  fossil  fuels,  less  pol¬ 
lution,  as  well  as  high  efficiency,  the  polymer  electrolyte  fuel  cell 
(PEFC)  is  regarded  as  one  of  the  most  promising  alternative  power 
sources  in  the  future.  It  is  expected  to  be  widely  employed  in  sta¬ 
tionary,  automotive  and  portable  sections.  However,  before  this 
can  occur,  several  technical  challenges  of  PEFCs  must  be  solved, 
such  as  cell  durability,  system  power  density,  fuel  storage,  genera¬ 
tion  and  delivery,  as  well  as  system  cost  to  ensure  a  proper  market 
penetration  [1-3]. 

In  a  single  PEFC  unit,  various  transport  processes  are  intricately 
coupled,  along  with  electrochemical  reactions  in  the  catalyst  layers. 
As  a  consequence,  water  and  heat  issues  are  always  ineluctable.  A 
typical  PEFC  consists  of  four  distinct  constituents,  namely,  bipo¬ 
lar  plate  (gas  channels  are  grooved  on  both  sides  of  bipolar 
plate),  gas  diffusion  layer  (including  the  micro  porous  layer),  cat¬ 
alyst  layer,  and  polymer  electrolyte  membrane.  On  one  hand,  the 
membrane  should  retain  high  water  content  to  transport  protons 
effectively  with  low  ohmic  resistance.  Hence,  gaseous  reactants 


Abbreviations:  GC,  gas  channel;  GDL,  gas  diffusion  layer;  CL,  catalyst  layer;  MEM, 
membrane;  PEFC,  polymer  electrolyte  fuel  cell. 
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(e.g.  H2,  02)  are  humidified  before  being  fed  into  fuel  cells.  On 
the  other  hand,  excessive  liquid  water  accumulation  within  fuel 
cells  would  block  reactant  pathways  to  reactive  sites  in  cata¬ 
lyst  layers,  resulting  in  the  so-called  flooding  situation.  Thus,  it 
is  evident  that  there  exist  two  conflicting  requirements  for  liquid 
water.  We  need  to  have  a  delicate  water  balance  inside  fuel  cells 
to  ensure  that  the  membrane  is  fully  hydrated  for  high  protonic 
conductivity,  while  severe  flooding  is  avoided,  especially  on  the 
cathode  side.  To  be  able  to  bring  about  this  balance,  a  profound 
understanding  of  water  transport  inside  fuel  cells  is  indispensable 
[4-7]. 

It  is  widely  recognized  that  the  flow  of  two  immiscible  phases 
(gas  and  liquid  water)  within  PEFCs  is  challenging.  While  PEFCs 
are  operating  under  certain  conditions  (e.g.  high  current  den¬ 
sities,  humid  environments,  and  cold  start-up),  liquid  water  is 
simultaneously  formed  in  all  components.  The  mechanisms  affect¬ 
ing  the  liquid  water  transport  are  distinct  in  different  layers.  A 
very  simple  categorization  of  the  two-phase  flow  in  PEFCs  can 
be  as  follows  [8,9]:  (1)  liquid  water  accumulation  and  transport 
in  the  CL,  (2)  two-phase  flow  in  the  GDL,  along  with  interfa¬ 
cial  coverage  at  the  GC-GDL  interface,  and  (3)  water  transport  in 
the  GC.  These  three  sub-processes  negatively  impact  the  perfor¬ 
mance  of  PEFCs.  For  instance,  in  the  CL,  excessive  liquid  water 
would  cover  active  catalyst  sites,  acting  as  an  additional  bar¬ 
rier  to  reactants  transport.  Based  on  the  preceding  descriptions 
of  water  transport  we  can  see  that  a  proper  water  management 
plays  a  central  role  in  the  development  and  commercialization  of 
PEFCs. 
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Nomenclature 

a  charge  transfer  coefficient;  net  water  transfer  coef¬ 

ficient 

Age  cross-sectional  area  of  gas  channel  (m2 ) 

Am  reactive  area  (m2) 

C  mass  fraction 

D  species  diffusivity  (m2  s-1 ) 

F  Faraday’s  constant,  96487  (C  mol-1 ) 

i0  exchange  current  density  (A  m-2 ) 

Jave  averaged  current  density  (A  m-2 ) 

j  volumetric  current  density  (A  m-3 ) 

/<o  intrinsic  permeability  ( m2 ) 

kr\  relative  permeability  for  liquid  phase 
krg  relative  permeability  for  gas  phase 

fccond  condensation  rate  (s-1 ) 

/<ev ap  evaporation  rate  (Pa-1  s-1 ) 

M  molecular  weight  (kg  mol-1) 

n  normal  direction  vector 

n\  exponent  for  the  effect  of  liquid  water  saturation  on 

species  diffusivity 

n2  exponent  for  relative  permeabilities 

n3  exponent  for  the  effect  of  liquid  water  saturation  on 

current  density 
p  pressure  (Pa) 

pg*,  water  vapor  pressure  (Pa) 

q  switch  function  for  phase  change  model 

R  universal  gas  constant  (8.134Jmol-1 1<-1);  mass 

source  due  to  phase  change  (kgm-3  s-1 ) 

RH  relative  humidity 

s  liquid  water  saturation 

5  source  term 

T  temperature  (I<) 

V  velocity  vector  (ms-1) 

U  intrinsic  velocity  vector  (ms-1) 

Y  molar  concentration  (mol  m-3 ) 

Greek  letters 

p  mass  density  (kg  m-3 ) 

s  porosity 

/x  dynamic  viscosity  (kg  m-1  s-1 ) 

r  stress  tensor  (N  m-2 ) 

o  surface  tension  (N  m-1 ) 

6  contact  angle  (°) 

a  cathode  catalyst  specific  area  (m2  m-3 ) 

pc  cathode  overpotential  (V) 

§  stoichiometric  ratio 

Superscripts  and  subscripts 
g  gas  phase 

1  liquid  phase 

m  momentum 

i  gaseous  species  index 

eff  effective  value 

c  capillary 

ref  reference  value 

im  immobile 

H20  water  vapor 

02  oxygen 

sat  saturation 

in  inlet 


Over  the  past  two  decades,  the  two-phase  flow  and  flooding 
phenomena  in  PEFCs  have  been  intensively  investigated  via  both 
experimental  [10-15]  and  numerical  methods  [16-25].  To  date, 
several  macroscopic  computational  fluid  dynamics  (CFD)  models 
for  the  two-phase  flow  in  PEFCs  are  available  in  literature,  which 
are  all  based  on  the  so-called  two-phase  Darcy’s  law  [26].  He  et  al. 
[27]  developed  a  two-phase  flow  model  for  the  cathode  GDL.  They 
solved  a  steady-state  transport  equation  of  liquid  water  flow  that 
derived  from  the  two-phase  Darcy’s  law.  The  equation  strongly 
resembled  a  general  scalar  transport  equation  with  convective  and 
diffusive  terms.  The  authors  also  assumed  both  capillary  diffusivity 
and  convective  coefficient  to  be  constants  for  numerical  stabil¬ 
ity.  Following  He’s  method,  Ye  and  Nguyen  [29]  also  derived  a 
similar  liquid  water  transport  equation.  Single-phase  flow  method¬ 
ology  was  employed  to  model  gas  flow,  and  two  phases  were 
coupled  by  phase  change.  The  effect  of  the  presence  of  liquid  water 
on  gas  flow  was  taken  into  account  only  by  correcting  gaseous 
species  diffusivities.  Another  popular  two-phase  model  for  PEFCs  is 
called  multiphase  mixture  (M2)  model,  which  has  been  employed 
widely  by  fuel  cell  researchers  [30-34].  Based  on  the  two-phase 
Darcy’s  law,  the  M2  model  for  multiphase,  multi-component  trans¬ 
port  in  capillary  porous  media  was  firstly  developed  by  Wang 
and  coworkers  [35,36].  However,  several  researchers  [28,37]  used 
volume-weighted  mixture  dynamic  viscosity  and  mass-weighted 
mixture  velocity  to  simplify  the  M2  model.  So,  the  applications  of 
this  modified  M2  model  would  be  limited.  Berning  et  al.  [38]  used 
the  so-called  multi-fluid  model  to  study  liquid  water  transport  in 
the  cathode  porous  layers.  Since  this  model  requires  a  multiphase 
solver,  and  needs  to  be  capable  of  coupling  species  transport,  phase 
change,  and  chemical  reactions  simultaneously,  it  entails  lots  of 
computational  efforts  and  is  prone  to  being  numerically  instable. 

In  an  operating  PEFC,  liquid  water  emerges  from  the  GDL  into 
the  GC,  in  the  form  of  small  droplets  and  slugs  [6].  These  droplets 
and  slugs  cover  the  GDL  surface  and  block  the  GC,  in  turn,  influence 
the  flooding  level  inside  the  diffusion  layers.  In  order  to  capture 
this  important  physical  phenomenon,  an  interactive  model  of  liq¬ 
uid  water  transport  between  the  GC  and  GDL  should  be  developed. 
In  most  previous  studies,  a  value  of  interfacial  saturation  or  cap¬ 
illary  pressure  at  the  GC-GDL  interface  was  specified.  Normally, 
this  value  was  assigned  to  zero,  corresponding  to  the  mist  flow 
assumption  in  the  GC.  However,  the  mist  flow  assumption  is  only 
valid  under  high  gas  flow  rates  in  the  GC,  which  are  not  encoun¬ 
tered  in  practice.  So  far,  only  a  few  researchers  have  numerically 
studied  the  water  coverage  effect  on  cell  performance.  Song  et  al. 
[39]  developed  a  one-dimensional  two-phase  analytical  model  to 
address  the  effect  of  liquid  water  saturation  at  the  GC-GDL  inter¬ 
face  on  the  transient  behavior  of  liquid  water  transport  inside  the 
cathode  GDL.  Results  showed  that  this  parameter  had  a  big  impact 
on  the  calculated  water  saturation  inside  the  GDL.  A  more  elaborate 
interfacial  coverage  model  was  proposed  by  Meng  and  Wang  [40]. 
In  their  work,  the  interfacial  liquid  saturation  at  the  GC-GDL  inter¬ 
face  was  assumed  to  be  a  simple  function  of  the  GDL  surface  contact 
angle,  current  density,  as  well  as  gas  inlet  velocity.  This  was  not 
based  on  a  derivation  and  they  also  used  the  mist  flow  assumption 
in  the  GC.  The  results  showed  that  the  interfacial  coverage  leaded  to 
higher  flooding  levels  inside  the  GDL  and  CL.  Berning  et  al.  [38]  used 
the  Hagen-Poiseuille  equation  to  relate  the  interfacial  water  satu¬ 
ration  with  local  pore  velocity  of  liquid  water  in  the  GDL.  Recently 
Basu  et  al.  [9]  proposed  to  apply  the  M2  model  into  the  GC  directly. 
In  this  approach,  the  GC  was  assumed  to  be  a  structured  porous 
medium,  and  then  the  two-phase  coupling  between  the  GDL  and 
GC  became  straightforward. 

In  this  work,  we  develop  a  two-phase  flow  model  for  the  cathode 
side  of  a  PEFC.  The  GC  is  assumed  to  be  a  structured  porous  medium 
with  the  porosity  of  1.0;  then,  the  two-phase  Darcy’s  law  is  applied 
to  the  GC.  Based  on  this  model,  we  study  the  liquid  water  flooding 
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Gas  flow  direction 
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Fig.  1.  (a)  Computational  domain  used  in  this  work;  (b)  and  the  corresponding  mesh  in  the  Y-Z  cross-section. 


in  the  GC  under  different  operating  conditions,  and  also  its  impact 
on  the  liquid  water  distribution  in  the  diffusion  layers.  Finally,  we 
address  the  effect  of  the  immobile  saturation  on  the  prediction  of 
the  liquid  water  distribution  in  the  diffusion  layers. 

2.  Two-phase  flow  model 

Fig.  1  shows  the  three-dimensional  computational  domain  and 
the  corresponding  mesh  in  the  Y-Z  cross-section.  In  this  work,  we 
focus  on  the  liquid  water  flooding  on  the  cathode  side  consisting  of 
the  CL,  GDL  and  GC.  For  the  sake  of  simplicity,  a  micro  porous  layer  is 
not  considered  here.  The  typical  cell  length  of  1 00  mm  is  used.  In  an 
operating  PEFC,  humidified  air  is  delivered  into  the  GC  and  diffuses 
through  the  GDL  into  the  CL.  There  oxygen  atoms  are  formed  on 
the  catalyst  where  they  combine  with  protons  from  the  membrane 
and  electrons  from  the  external  circuit  to  generate  liquid  water  and 
heat.  The  following  main  assumptions  are  made  in  our  model  for 
numerical  simplicity: 

(1)  Ideal  gas  mixture 

(2)  Isotropic  and  homogeneous  diffusion  layers 

(3)  Laminar  flow  due  to  small  Reynolds  numbers  and  pressure  gra¬ 
dients 

(4)  Isothermal  condition. 

It  is  well  known  and  proofed  experimentally  that  the  GDL  is 
anisotropic  when  comparing  in-plane  and  through-plane  proper¬ 
ties  like  thermal  and  electrical  conductivities,  species  diffusivities 
and  permeabilities  [41  ].  Anisotropic  thermal  and  electrical  conduc¬ 
tivities  do  not  play  a  role  in  our  model,  since  heat  and  electron 
transport  is  excluded.  According  to  Tomadakis  and  Sotirchos  [42], 
the  in-plane  diffusivity  is  around  20%  larger  than  the  through-plane 
one.  The  conclusion  indicates  that  the  disregard  of  the  anisotropic 
diffusivity  in  this  study  should  not  impact  the  general  tendencies  of 
the  results.  The  same  argument  holds  for  the  permeabilities,  since 
the  through-plane  permeability  of  carbon-based  GDL  is  around  50% 
bigger  than  the  in-plane  one  [43]. 

2.1.  Governing  equations  of  gas  phase 

The  mass  conservation  equation  of  the  gas  mixture  is  given  as 
follows: 

V.(pgVg)  =  Sg  (1) 

where  Vg  denotes  the  superficial  or  Darcy  velocity  of  the  gas  phase, 
pg  is  the  density  of  the  gas  mixture,  and  Sg  is  the  mass  source  term 
due  to  the  chemical  reaction  and  phase  change. 


The  Darcy  velocity  is  related  to  the  intrinsic  or  pore-scale  veloc¬ 
ity  by  the  following  definition: 

Vg  =  (l  -s)eUg  (2) 

where  Ug  denotes  the  intrinsic  velocity  of  the  gas  phase,  s  is  the 
porosity,  and  s  represents  the  liquid  water  saturation  defined  as 
the  volume  faction  of  the  pores  occupied  by  the  liquid  water. 

The  steady-state  momentum  equation  of  the  gas  phase  can  be 
derived  based  on  the  two-phase  flow  methodology  [44],  which  is 
given  as: 

-J-L-?v.(„sv>i;i).-yPg+slLJ5v.(T)+s™  (3) 

where  pg  is  the  gas  phase  pressure,  r  is  the  stress  tensor,  and  Sm  rep¬ 
resents  the  Darcy  source  term  accounting  for  the  viscous  resistance 
imposed  by  the  pore  structure  of  porous  media  and  the  presence 
of  liquid  water  in  the  pore  space. 

The  above-listed  mass  and  momentum  conservation  equations 
of  the  gas  phase  are  applied  to  the  whole  computational  domain. 
Since  we  regard  the  GC  as  a  structured  porous  medium  in  this  work, 
the  porosity  of  the  GC  is  assumed  to  be  unity.  It  is  worth  noting 
that  the  two-phase  Darcy’s  law  for  the  gas  phase  can  be  obtained 
by  neglecting  the  inertia  and  diffusive  terms  in  Eq.  (3).  Flowever,  to 
keep  the  same  order  of  the  momentum  equations  used  in  the  GC 
and  diffusion  layers,  Eq.  (3)  is  adopted  in  this  work. 

In  automotive  applications,  we  have  water  vapor,  oxygen,  and 
nitrogen  species  on  the  cathode  side.  The  species  transport  equa¬ 
tion  can  be  expressed  as: 

V  .  (pgQVg)  =  V  .  (pgDfVCi)  +  Si  (4) 

in  which  Q  denotes  the  mass  fraction  of  each  component,  is  the 
effective  diffusivity  accounting  for  the  presence  of  liquid  water  and 
the  pore  structure  of  porous  media,  and  S2-  represents  the  species 
source  term  owing  to  the  oxygen  reduction  reaction  (ORR)  and 
phase  change.  In  our  calculations,  nitrogen  is  regarded  as  an  inert 
component. 

2.2.  Liquid  water  transport  equation 

Liquid  water  is  assumed  to  be  a  single-component  phase,  and 
the  gaseous  species  diffusion  in  the  liquid  water  is  neglected.  This 
is  a  valid  assumption  as  a  first  approach  since  the  liquid  diffusivities 
are  several  orders  of  magnitudes  lower  than  the  gaseous  diffusiv¬ 
ities.  The  mass  conservation  equation  of  the  liquid  water  is  given 
as: 

V-(AVJ)  =  Si  (5) 
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where  p\  is  the  liquid  water  density,  5!  is  the  mass  source  term, 
and  V\  denotes  the  Darcy  velocity  of  the  liquid  water,  which  can  be 
given  by  the  two-phase  Darcy’s  law: 

Vi  =  -^VPl  (6) 

fi\ 

Recalling  Eq.  (3),  we  can  safely  drop  the  inertia  and  diffusive 
terms  due  to  their  extremely  small  magnitudes  compared  to  the 
Darcy  term.  Then  the  two-phase  Darcy  equation  for  the  gas  phase 
is  obtained  as: 

vg  =  -^vpg  (7) 

In  Eqs.  (6)  and  (7),  p\  denotes  the  liquid  water  pressure,  k0  denotes 
the  intrinsic  permeability  of  porous  media,  /i\  and  /xg  are  the 
dynamic  viscosities  of  the  liquid  water  and  gas,  respectively,  kT\ 
represents  the  relative  permeability  for  the  liquid  phase,  and  krg 
represents  the  relative  permeability  for  the  gas  phase.  The  primary 
difference  between  these  two  equations  and  the  single-phase  form 
stems  from  the  presence  of  the  relative  permeability,  which  desig¬ 
nates  the  influence  of  the  reduced  void  space  for  each  phase  due 
to  the  existence  of  the  other  phase.  Finally,  the  macroscopic  cap¬ 
illary  pressure  is  introduced  to  represent  the  pressure  difference 
between  the  gas  and  liquid  phases,  given  as: 


Pc(s)  =  Pg-Pi  (8) 

Eqs.  (5)-(8)  can  be  combined,  and  after  a  few  algebraic  manipula¬ 
tions  the  liquid  water  transport  equation  is  obtained  as: 


Here,  the  gas  drag  effect  (the  term  on  the  left-hand  side)  on  the 
liquid  water  movement  and  the  capillary  diffusion  (first  term  on 
the  right-hand  side)  are  both  included. 


2.3.  Constitutive  correlations  and  source  terms 


To  close  the  governing  equations,  several  constitutive  corre¬ 
lations  are  needed,  such  as  state  equation,  capillary  pressure- 
saturation  and  relative  permeability-saturation  relationships.  Ideal 
gas  mixture  law  is  used  to  calculate  the  composition-dependent  gas 
phase  density  as  follows: 


Pg  - 


Pg 

RTZfi/Mt 


(10) 


in  which,  pg  refers  to  the  gas  pressure,  Q  is  the  mass  fraction  of  each 
component,  and  R  is  the  universal  gas  constant.  Due  to  the  small 
pressure  gradient,  the  mass  density  of  the  liquid  water  is  assumed 
to  be  constant. 

The  multi-component  diffusion  in  the  gas  phase  can  be 
described  by  Stefan-Maxwell  equation.  For  simplicity,  the  Fields 
law  is  commonly  used  [45].  Taking  into  account  the  presence  of 
liquid  water  and  the  pore  structure  of  porous  media,  the  effective 
species  diffusivity  can  be  given  based  on  Bruggeman  correlation: 


Dfff  =  (1  -spS13Di 


(11) 


where  D*  is  the  molecular  diffusivity  that  depends  on  the  tem¬ 
perature  and  pressure.  For  water  vapor  and  oxygen,  we  have  the 
following  empirical  expressions  [29]: 


(12) 

(13) 


Regarding  the  capillary  pressure,  we  know  that  it  is  depen¬ 
dent  on  the  interface  curvature  at  the  micro  scale.  We  traditionally 
presume  that  capillary  pressure  is  solely  a  function  of  saturation 
at  the  macro  scale,  as  expressed  in  Eq.  (8).  However,  the  hys¬ 
teresis  phenomenon  of  capillary  pressure-saturation  curves  during 
consecutive  drainage  and  imbibition  processes  in  porous  media 
indicates  that  other  variables  may  also  influence  the  capillary  pres¬ 
sure.  Hassanizadeh  and  Gray  [46]  claim  that,  based  on  a  rigorous 
thermo-dynamical  derivation  the  capillary  pressure  is  not  only  a 
function  of  saturation  but  also  of  the  specific  areas  of  the  three 
interfaces.  However,  to  date,  this  model  is  still  hard  to  be  inte¬ 
grated  with  current  CFD  codes,  since  it  involves  several  additional 
variables.  There  is  also  another  weakness  in  the  traditional  capil¬ 
lary  pressure-saturation  relationship,  which  ignores  the  dynamic 
effect  [47]  in  the  pressure  difference  between  two  phases  under 
unsteady  state  situations.  In  this  work,  we  still  follow  the  traditional 
approach,  and  use  the  well-known  Leverett  function,  expressed  as: 


/  £  \  1/2 

Pc  =a  cos  0(j^)  J(s) 

f  1.417(1  —  s)  —  2.120(1  —  s)2  + 1.263(1  -s)3  for  0  <  90° 
[^  1.417s  -  2.120s2  +  1.263s3  for  0  >  90° 


(14) 

(15) 


As  to  the  relative  permeabilities,  power-form  correlations  are 
commonly  used  in  the  absence  of  experimental  data.  They  orig¬ 
inally  come  from  sand/rock-type  porous  media  with  a  typical 
porosity  of  0.1 -0.4  [29].  Therefore,  the  relative  permeabilities  for 
both  gas  and  liquid  phases  can  be  expressed  as: 


krg  —  (1  —  S)  2 
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(17) 

(18) 


In  these  equations,  n2  is  the  material  coefficient,  seff  is  called 
the  effective  saturation,  and  sim  represents  the  immobile  water 
saturation. 

The  liquid  water  flux  is  assumed  to  be  zero,  before  it  forms  con¬ 
ducting  pathways  in  the  porous  media  of  interest.  As  indicated  in 
Eq.  (17),  when  the  liquid  water  saturation  is  less  than  the  thresh¬ 
old  value  sim,  the  relative  permeability  for  the  liquid  phase  is  zero; 
thus,  we  obtain  the  zero  water  flux.  Normally,  the  immobile  sat¬ 
uration  not  only  depends  on  the  material  properties  of  porous 
media,  but  also  can  be  affected  by  the  flow  conditions  (such  as 
phase  change,  boundary  condition,  and  flow  history).  Only  a  few 
researchers  [38,48]  have  investigated  the  impact  of  this  parame¬ 
ter  on  the  two-phase  flow  in  the  diffusion  layers  of  a  PEFC  using 
one-dimensional  simulations. 

The  simplified  Tafel  equation  is  employed  to  describe  the  rela¬ 
tively  sluggish  ORR  in  the  cathode  CL  [22]: 


jc  =  (1 


-s)"3aiQef 
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^02 .  ref 


fig 
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(19) 


Here,  the  liquid  water  coverage  effect  is  considered  by  the 
exponent  n3,  Y02,ref  denotes  the  reference  molar  concentration  of 
oxygen,  and  r/c  is  the  cathode  overpotential. 

The  phase  change  between  water  vapor  and  liquid  water  is  con¬ 
sidered  by  a  nonequilibrium  phase  change  model  [49],  which  is 
expressed  as: 


Ki  =  K 
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i  +  |yH2oPg  -  ph2o  |  /(^oPg  -  p^o)  , 

In  Eq.  (20),  pj^  denotes  the  water  vapor  pressure,  which  can  be 
obtained  from  the  following  empirical  formula  as  a  function  of 
temperature  [44]: 

p??'  /101,325 

logjQ  0  =  -2.1794  +  0.02953(7  -  273.15) 

-9.1837  x  10_5(T  —  273. 15)2 
+  1.4454  x  10“7(T-273.15)3  (22) 

The  source  terms  of  the  conservation  equations  in  different 
regions  on  the  cathode  side  are  listed  in  Table  1.  Since  the  mem¬ 
brane  is  not  included,  the  liquid  water  exchange  between  the  CL 
and  membrane  is  accounted  for  by  means  of  the  net  water  trans¬ 
fer  coefficient.  Note  that  the  generated  water  in  the  CL  is  assumed 
to  be  in  the  liquid  form  due  to  the  nanostructure  of  the  CL.  In  this 
work,  we  assume  the  GC  to  be  a  structured  porous  medium  asso¬ 
ciated  with  a  given  intrinsic  permeability,  which  can  be  obtained 
from  numerical  experiments  or  Hagen-Poiseuille  law  [50]. 


Table  2 

Geometrical,  physical  and  operating  parameters. 


Parameter 

Value 

GC  width/height/length 

1.0/0.5/100  mm 

Shoulder  width 

1.0  mm 

GDL/CL  thickness 

0.15/0.015  mm 

GC/GDL/CL  porosity,  e 

1.0/0.6/0.4 

GC/GDL/CL  intrinsic  permeability,  ko 

1.4  x  10_8/3.  x  10-12/3.  x  10“14 

Cell  temperature,  T 

353.15  K 

Operating  pressure,  pref 

2  atm 

n-i/n3 

2.0  /2.0  [9] 

GDL  /CL  /GC  n2 

4.0/4.0/5.0  [9] 

Liquid  water  density,  p\ 

972  kg  m-3 

Dynamic  viscosity  of  liquid  water/gas 

3.5  x  10“4/2.03  x  10“5  Pas  [29] 

mixture,  ii 

Evaporation  rate,  kev ap 

5.0  x  10-5  Pa-1  s-1  [29] 

Condensation  rate,  /<cond 

1.0s-1  [29] 

GC/GDL/CL  contact  angle,  6 

60/120/110° 

Surface  tension,  liquid-water-air,  cr 

0.0625  Nm-1 

Charge  transfer  coefficient,  ac 

1.0 

Reference  current  density  x  ratio  of 

2.0  x  104  Am-3 

reaction  surface  to  catalyst  volume 
on  cathode  side,  aiTf[ 

Reference  oxygen  molar  concentration, 

40.88  mol  m-3 

Yo’ 

Net  water  transfer  coefficient,  a 

0.0 

3.  Results  and  discussion 


2.4.  Boundary  conditions  and  numerical  implementation 


The  gas  inlet  velocity  on  the  cathode  side  can  be  determined 
by  the  specified  stoichiometric  ratio  and  the  calculated  average 
current  density  /aVe*  as  follows: 


(23) 


Geometrical,  physical,  and  operating  parameters  are  summa¬ 
rized  in  Table  2.  We  assume  uniform  distribution  of  the  cathode 
overpotential  in  the  CL,  and  the  average  current  density  can  be 
computed  from  Eq.  (23).  In  what  follows,  we  highlight  the  impor¬ 
tance  of  modeling  GC  flooding  in  the  numerical  studies  of  PEFCs. 
And,  its  impact  on  the  cell  performance  and  liquid  water  distribu¬ 
tion  is  explored  in  detail.  At  last,  we  address  the  effect  of  immobile 
saturation  on  the  liquid  water  distribution  in  diffusion  layers. 
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(24)  3.1.  Modeling  of  GC  flooding 


where  Am  and  Ac  are  the  reactive  and  GC  cross-sectional  areas, 
respectively,  and  Yin02  denotes  the  inlet  molar  concentration  of 
oxygen,  which  is  determined  by  the  inlet  relative  humidity  RH  and 
local  gas  pressure,  expressed  as: 


RH-pffpN 
Pg  ) 


(25) 


In  experiments  with  PEFCs,  a  large  amount  of  water  droplets  and 
slugs  forming  in  the  GC  are  observed,  which  are  finally  removed 
out  of  the  channel  mainly  due  to  the  gas  drag  force.  In  previous 
numerical  studies,  two-phase  flow  in  the  GC  was  overwhelmingly 
simplified  by  the  mist  flow  assumption.  In  this  work,  we  apply  the 
two-phase  Darcy’s  law  to  the  GC,  and  investigate  its  effect  on  the 
cell  performance  and  liquid  water  distribution.  Fig.  2  shows  two 
polarization  curves  (cathode  potential  vs.  current  density)  under 


The  gas  pressure  at  the  GC  outlet  is  assumed  to  be  2  atm.  The 
symmetric  boundary  conditions  are  employed  at  the  sidewalls  of 
the  GDL  and  CL  due  to  the  repeated  structure  in  a  PEFC  stack.  The 
slip  and  impermeable  wall  conditions  are  specified  at  the  GC  walls, 
since  the  GC  is  assumed  to  be  a  structured  porous  medium.  The 
wall  resistance  is  lumped  into  the  Darcy’s  source  term  as  shown 
in  Table  1.  For  the  remaining  boundaries,  the  conditions  of  no-slip 
and  impermeable  wall  are  imposed. 

The  set  of  steady-state  governing  equations  and  boundary  con¬ 
ditions  given  above  are  discretized  using  finite  volume  method 
with  second-order  schemes  based  on  the  commercial  CFD  solver 
FLUENT.  With  the  aid  of  user-defined  functions,  the  liquid  water 
transport  equation,  source  terms  and  constitutive  correlations 
are  coupled  with  the  general  multi-component  gas  transport 
equations.  The  SIMPLE  (semi-implicit  method  for  pressure-linked 
equations)  algorithm  is  utilized  to  couple  pressure  and  velocity,  and 
the  AMG  (algebraic  multi-grid)  method  in  conjunction  with  Gauss- 
Seidel  type  smoother  is  used  to  solve  resultant  nonlinear  algebraic 
equations.  In  all  simulations  presented  in  the  following  section,  the 
scaled  values  of  the  equation  residuals  are  smaller  than  10-6. 


Fig.  2.  Polarization  curves  (cathode  potential  vs.  current  density)  taking  into 
account  GC  flooding,  and  neglecting  GC  flooding.  It  is  clearly  seen  that  the  cell 
performance  reduces  due  to  the  liquid  water  flooding  in  the  GC. 
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Table  1 

Source  terms  of  the  conservation  equations  in  different  regions  on  the  cathode  side. 
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the  situations  of  considering  the  GC  flooding,  and  neglecting  the 
GC  flooding,  separately.  It  can  be  seen  that  neglecting  the  GC  flood¬ 
ing  would  overestimate  the  cell  performance.  We  also  validate  our 
model  against  experimental  data  in  terms  of  gas  pressure  drop 
in  the  GC  [6,9].  Based  on  an  operating  PEFC,  Hussaini  et  al.  [6] 
measured  the  gas  pressure  drops  under  different  flooding  situa¬ 
tions  in  the  GC.  In  our  simulations,  the  same  GC  dimensions  (see 
Table  2)  and  operating  conditions  (RH  =  66%)  as  in  Ref.  [6]  are  used. 
Comparisons  of  numerical  and  experimental  results  are  shown 
in  Fig.  3.  The  excellent  match  is  obtained  at  the  current  density 
of  0.5  A  cm-2.  Increasing  the  stoichiometric  ratio  does  not  invoke 
obvious  increase  of  the  gas  pressure  drop,  since  more  liquid  water 
can  be  swept  out  of  the  channel  under  the  higher  gas  flow  rate. 
However,  our  model  overpredicts  the  gas  pressure  drops  at  the 
current  density  of  0.8  A  cm-2,  and  underestimates  the  gas  pressure 
drops  at  the  low  current  density  of  0.2  A  cm-2,  respectively.  This 
may  be  partly  due  to  the  unresolved  water  exchange  between  the 
cathode  and  anode  sides  in  our  half  cell  model.  It  also  needs  to 
point  out  that  our  channel  flooding  model  is  capable  of  predicting 
the  tendency  of  gas  pressure  drop  as  increasing  the  stoichiometric 
ratio,  as  shown  in  Fig.  3. 

Fig.  4  shows  the  comparison  of  the  liquid  water  distributions  at 
the  middle  cross-section  (Y=  1  x  10-3  m)  of  the  cathode  side  under 
two  situations:  considering  the  GC  flooding  (Fig.  4a),  and  neglecting 
the  GC  flooding  (Fig.  4b).  The  corresponding  operating  conditions 
are  as  follows :  RH  =  66%,  £c  =  2.0,  sim  =  0,  and  JaVe  =  0.2  A  cm-2.  Much 
more  liquid  water  is  found  in  the  diffusion  layers  due  to  the  water 
coverage  effect  at  the  GC-GDL  interface.  The  cathode  overpotential 
increases  to  0.392  V  when  taking  into  account  the  GC  flooding.  This 
indicates  the  decreased  cell  performance. 

At  the  same  operating  conditions,  the  liquid  water  distributions 
at  the  middle  cross-section  (Z  =9.0x1 0-4  m)  of  the  GDL  are  plotted 
in  Fig.  5.  The  liquid  water  saturation  at  the  inlet  portions  is  zero 
due  to  the  partially  humidified  air  feeding.  When  neglecting  the  GC 
flooding  (Fig.  5b)  more  liquid  water  is  found  under  the  land,  this  is 
because  the  liquid  water  under  the  land  has  to  travel  a  somewhat 


Stoichiometric  ratio,  & 

Fig.  3.  Gas  pressure  drop  validation  at  different  current  densities  and  stoichiometric 
ratios  (experimental  data  are  extracted  from  Ref.  [9]). 
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Fig.  4.  Distributions  of  the  liquid  water  saturation  at  the  middle  cross-section 
( Y=  1 .0  x  1 0-3  m)  of  the  cathode  side  (operating  conditions:  RH  =  66%,  £c  =  2.0,  Sim  =  0, 
and  /ave  =  0.2  A  cm-2 ). 


longer  path  to  reach  the  GC.  The  liquid  water  saturation  reduces 
slightly  along  the  flow  direction  owing  to  the  decreased  current 
density.  When  the  GC  flooding  model  is  included,  the  liquid  water 
saturation  increases  along  the  flow  direction  (Fig.  5a).  This  opposite 
water  distribution  indicates  that  the  GC  flooding  impacts  the  liquid 
water  distribution  in  the  diffusion  layers  considerably.  Neglecting 
the  GC  flooding  would  lead  to  the  incorrect  prediction  of  liquid 
water  distribution  in  the  diffusion  layers.  In  addition,  the  liquid 
water  distribution  in  the  transverse  direction  of  the  GC  is  more 
uniform,  which  is  also  caused  by  the  presence  of  liquid  water  in 
the  GC. 

Fig.  6a  shows  the  liquid  water  distributions  at  the  middle  cross- 
section  (X=  0.05  m)  of  the  cathode  side.  Liquid  water  flooding  in 
the  GC  results  in  almost  even  distribution  of  the  liquid  water  at  the 
cross-section  (left).  However,  a  distinct  gradient  of  the  liquid  water 
saturation  is  found  when  we  neglect  the  liquid  water  flooding  in  the 
GC  (right).  Fig.  6b  displays  the  distributions  of  the  mole  fraction 
of  oxygen  at  the  same  cross-section.  Lower  oxygen  concentrations 
are  found  under  the  land  due  to  the  transport  limitation.  When 
considering  the  GC  flooding,  much  more  liquid  water  is  present  in 


With  GC  flooding  model 


Fig.  5.  Distributions  of  the  liquid  water  saturation  at  the  middle  cross-section 
(Z=9.0  x  10_5m)  of  the  GDL  (operating  conditions:  RH  =  66%,  £c  =  2.0,  sim  =  0,  and 
Lve  =  0.2  A  cm-2 ). 
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Without  GC  flooding  model  Without  GC  flooding  model 


(a) 


(b) 


Fig.  6.  Distributions  of  (a)  the  liquid  water  saturation,  (b)  the  mole  fraction  of  oxygen 
at  the  middle  cross-section  (X=0.05  m)  of  the  cathode  side  (operating  conditions: 
RH  =  66%,  £c  =  2.0,  sim  =  0,  and  Jave  =  0.2  A  cm-2 ). 


the  diffusion  layers  (Fig.  6a).  This  gives  rise  to  the  higher  oxygen 
diffusion  resistance.  Therefore,  even  lower  oxygen  concentration 
is  observed  (Fig.  6b,  left  one). 

Fig.  7  illustrates  the  influence  of  the  stoichiometric  ratio 
on  the  liquid  water  distribution  at  the  middle  cross-section 
(Z=4.15  x  10-4  m)  of  the  GC.  The  operating  conditions  are: 
RFI  =  66%,  sim  =  0,  and  JaVe  =  0.5  A  cm-2.  When  the  applied  stoichio¬ 
metric  ratio  equals  to  2.0,  the  highest  liquid  water  saturation  is 
predicted  to  be  0.23  at  the  end  of  the  GC  (Fig.  7a).  We  also  can 
observe  the  water  vapor  front  at  which  the  water  vapor  reaches  its 
saturated  value.  With  the  increase  of  the  stoichiometric  ratio,  the 
liquid  water  flooding  is  mitigated  in  the  GC.  This  can  be  explained 
by  the  fact  that  more  liquid  water  is  removed  out  of  the  GC  with 
the  increase  of  the  gas  flow  rate.  Meanwhile,  the  water  vapor  front 
moves  fast  towards  the  end  of  the  GC,  which  indicates  that  more 
water  is  transported  out  of  the  GC  in  the  water  vapor  form.  When 
the  applied  stoichiometric  ratio  equals  to  4.0,  only  a  little  liquid 
water  is  left  at  the  end  of  the  GC.  And,  the  cathode  overpotential 
decreases  to  0.348  V,  indicative  of  the  increased  cell  performance. 

Fig.  8  illustrates  the  influence  of  the  inlet  air  relative  humid¬ 
ity  on  the  liquid  water  distribution  at  the  middle  cross-section 
(Z=4.15  x  10-4  m)  of  the  GC.  Three  different  values  of  relative 


Fig.  7.  Distributions  of  the  liquid  water  saturation  at  the  middle  cross-section 
(Z  =  4.15  xlO-4)  of  the  GC  at  different  stoichiometric  ratios:  (a)  £c  =  2.0,  (b)  £c  =  3.0, 
(c)  and  £c  =  4.0  (operating  conditions:  RH  =  66%,  sim  =  0,  and  /ave  =  0.5  A  cm-2 ). 
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Fig.  8.  Distributions  of  the  liquid  water  saturation  at  the  middle  cross-section 
(Z  =  4.15  x  10-4)  of  the  GC  at  different  inlet  relative  humidity:  (a)  RH  =  46%, 
(b)  RH  =  66%,  (c)  and  RH  =  86%  (operating  conditions:  £c  =  2.0,  sim  =  0,  and 
/ave  =  0.8  A  cm-2). 


humidity  are  assigned,  namely,  46%,  66%,  and  86%.  The  operating 
conditions  are:  £c  =  2.0,  sim  =  0,  and  JaVe  =  0.8  A  cm-2.  At  the  low  rel¬ 
ative  humidity  of  46%  (Fig.  8a),  only  less  than  half  GC  is  flooded 
by  the  liquid  water,  and  most  generated  water  is  removed  out  of 
the  GC  in  the  water  vapor  form  as  mentioned  above.  As  a  conse¬ 
quence,  we  get  the  minimum  cathode  overpotential  of  0.376  V  in 
all  three  cases.  However,  it  is  worth  noting  that  low  relative  humid¬ 
ity  would  result  in  the  membrane  dehydration  mainly  at  the  inlet 
portion,  which  decreases  the  cell  performance.  From  Fig.  8,  we  also 
can  see  that  more  and  more  liquid  water  accumulates  in  the  GC  with 
the  increase  of  the  inlet  relative  humidity.  Although  increasing  the 
inlet  relative  humidity  can  slightly  increase  the  gas  flow  rate  at  the 
same  current  density  and  stoichiometric  ratio,  the  phase  change 
dominates  the  liquid  water  flooding  in  the  GC.  Therefore,  we  obtain 
the  maximum  water  saturation  up  to  0.24,  when  the  inlet  relative 
humidity  equals  to  86%  (Fig.  8c).  And,  almost  all  the  GC  is  flooded 
by  the  liquid  water.  Comparison  between  Fig.  8b  and  Fig.  7a  shows 
that  increasing  the  current  density  from  0.5  to  0.8  A  cm-2  gives  the 
indiscernible  change  of  the  GC  flooding.  This  is  attributed  to  the 
fact  that  the  high  current  density  corresponds  to  the  high  gas  flow 
rate,  which  can  drag  more  liquid  water  out  of  the  channel. 

In  this  work,  we  apply  Eq.  (9)  to  track  the  liquid  water  flooding 
in  the  GC  as  a  first  attempt.  By  examining  Eq.  (9),  we  can  see  that 
two  mechanisms  determine  the  liquid  water  transport  in  the  GC, 
namely,  gas  drag  force  and  capillary  action.  However,  their  magni¬ 
tudes  strongly  depend  on  the  assumed  material  properties,  such  as 
GC  contact  angle  and  relative  permeability  for  each  phase.  In  what 
follows,  we  evaluate  the  sensitivity  of  the  simulated  GC  flooding  to 
these  assumed  material  properties.  Fig.  9  shows  the  contact  angle 
effect  on  the  liquid  water  distribution  at  the  middle  cross-section 
(Z=4.15  x  10-4  m)  of  the  GC.  It  needs  to  note  that  the  real  GC  nor¬ 
mally  comprises  of  the  hydrophilic  sidewalls  and  is  enclosed  by 
the  hydrophobic  GDL  surface.  However,  in  this  work,  we  lump  all 
the  wettability  features  to  a  GC  contact  angle.  From  Fig.  9  it  can 
be  seen  that  the  liquid  water  flooding  in  the  GC  is  insensitive  to 
the  contact  angle  at  the  fixed  cathode  overpotential  of  0.38  V.  This 
indicates  that  the  capillary  action  is  not  the  dominant  mechanism 
of  the  liquid  water  transport  in  the  GC. 

Fig.  10  demonstrates  the  sensitivity  of  the  liquid  water  flooding 
to  the  gas  drag  coefficient  in  the  GC.  According  to  Eq.  (9),  we  can 
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Fig.  9.  Distributions  of  the  liquid  water  saturation  at  the  middle  cross-section 
(Z=4.15  x  10“4)  of  the  GC  at  different  GC  contact  angles:  (a)  0  =  60°,  (b)  0  =  80°,  and 
(c)  0  =  110°  (operating  conditions:  RH  =  66%,  £c  =  2.0,  sim  =  0,  and  r]c  =  0.38  V). 


see  that  the  gas  drag  effect  is  directly  determined  by  the  relative 
permeabilities  for  both  phases  in  the  GC.  By  adjusting  the  exponent 
n2  in  the  definition  of  the  relative  permeability,  we  can  evaluate 
the  gas  drag  effect  on  the  liquid  water  flooding  in  the  GC.  Fig.  10 
shows  that  the  GC  flooding  is  very  sensitive  to  the  applied  material 
property  n2.  With  the  decrease  of  n2,  more  and  more  liquid  water  is 
swept  out  of  the  channel  by  the  gas  flow.  This  indicates  that  proper 
estimation  of  the  gas  drag  force  is  crucial  to  the  GC  flooding  model 
used  in  this  work.  When  the  GC  is  assumed  to  be  a  structured  porous 
media,  the  nominal  relative  permeability  for  each  phase  cannot  be 
only  a  function  of  water  saturation,  but  also  of  other  variables  (e. 
g.  gas  flow  rate  and  current  density).  Maybe  experimental  data  can 
help  fit  an  empirical  formulation. 

3.2.  Immobile  saturation  effect 

The  immobile  saturation  could  be  an  important  parameter  in  the 
modeling  of  liquid  water  transport  in  the  diffusion  layers.  Before  the 
liquid  water  forms  conducting  pathways  (i.e.  the  water  saturation 
is  less  than  sim),  its  relative  permeability  equals  to  zero;  otherwise, 


n2  =  5.0  (for  GC)  1^  =  0.479  A  cm  2 


Flow  direction 


Fig.  10.  Impact  of  the  gas  drag  coefficient  on  the  water  saturation  distributions  at 
the  middle  cross-section  (Y=  1.0  x  10-3  m)ofthe  cathode  side  (operating  conditions: 
RH  =  66%,  £c  =  2.0,  sim  =  0,  and  r)c  =  0.36  V). 


Liquid  water  saturation  increases 


Fig.  11.  Effect  of  the  immobile  saturation  on  the  liquid  water  distribution  at  the 
middle  cross-section  (Z=  9.0  x  10-5  m)  of  the  GDL  (operating  conditions:  RH  =  66%, 
£c  =  2.0,  and  r]c  =  0.36  V). 


the  effective  saturation  is  employed  in  the  relative  permeability- 
saturation  correlation.  Since  the  gas  phase  is  always  continuous 
throughout  the  whole  domain,  its  relative  permeability  is  given 
by  Eq.  (16).  Ju  [48]  used  the  effective  saturation  in  the  capillary 
pressure-saturation  correlation.  However,  it  is  better  to  work  with 
the  actual  water  saturation  s  as  shown  in  Eq.  (15). 

As  stated  earlier,  the  immobile  saturation  depends  on  the 
flow  conditions  and  microstructure  of  porous  media,  normally  an 
approximate  value  is  assigned  to  it.  In  this  work,  we  assign  a  con¬ 
stant  value  of  0.05  to  the  immobile  saturation,  and  investigate  its 
effect  on  the  liquid  water  distribution  in  the  diffusion  layers.  It 
needs  to  point  out  that  when  the  immobile  saturation  is  taken  into 
account,  the  exponent  n2  should  be  corrected  to  fit  experimental 
data.  However,  due  to  the  absence  of  experimental  data,  we  keep 
it  unchanged  as  a  first  attempt. 

Fig.  1 1  shows  the  effect  of  the  immobile  saturation  on  the  liq¬ 
uid  water  distribution  at  the  middle  cross-section  (Z=9.0  x  10-5) 
of  the  GDL.  The  operating  conditions  are:  RH  =  66%,  £c  =  2.0,  and 
r] c  =  0.36  V.  When  considering  the  immobile  saturation  (lower  one), 
more  liquid  water  resides  in  the  diffusion  layers.  As  a  result,  the 
current  density  decreases  slightly.  Unlike  the  GC  flooding,  the 
immobile  saturation  does  not  change  the  tendency  of  the  liquid 
water  distribution  along  the  flow  direction.  It  only  increases  the 
flooding  level  in  the  diffusion  layers. 


4.  Conclusions 

A  two-phase  flow  model  for  the  cathode  side  of  a  PEFC  has  been 
developed  in  this  work.  The  salient  feature  is  that  we  assume  the 
GC  to  be  a  structured  porous  medium;  then,  the  liquid  water  flood¬ 
ing  in  the  GC  can  be  modeled  by  the  two-phase  Darcy’s  law  as  a 
first  attempt.  We  investigate  the  GC  flooding  under  different  oper¬ 
ating  conditions,  and  its  impact  on  the  liquid  water  distribution  in 
the  diffusion  layers.  We  also  study  the  effect  of  the  immobile  sat¬ 
uration  on  the  predicted  liquid  water  distribution  in  the  diffusion 
layers.  The  main  findings  and  conclusions  are  summarized  in  the 
following: 

(1)  Neglecting  the  GC  flooding  leads  to  the  incorrect  prediction  of 
the  liquid  water  distribution  in  the  diffusion  layers,  and  also 
overestimates  the  cell  performance. 

(2)  Increasing  the  stoichiometric  ratio  does  not  invoke  obvious 
increase  of  the  gas  pressure  drop  in  the  GC,  since  more  liquid 
water  can  be  removed  out  of  the  channel  under  the  higher  gas 
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flow  rate.  This  indicates  that  practically  we  can  optimize  the 
gas  flow  rate  to  achieve  the  best  cell  performance. 

(3)  Decreasing  the  relative  humidity  of  the  inlet  gas  can  mitigate 
the  GC  flooding,  since  more  liquid  is  transported  out  of  the 
channel  in  the  water  vapor  form.  However,  a  proper  selection 
of  relative  humidity  is  needed  to  balance  two  requirements: 
preventing  the  membrane  dehydration,  and  mitigating  the  GC 
flooding. 

(4)  The  gas  flow  is  the  main  driving  force  for  the  liquid  water  trans¬ 
port  in  the  GC.  The  liquid  water  flooding  in  the  GC  is  insensitive 
to  the  capillary  action. 

(5)  When  considering  the  immobile  saturation  in  the  model,  more 
liquid  water  is  predicted  to  be  in  the  diffusion  layers,  and  the 
cell  performance  decreases  slightly. 
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